Datos abiertos SSA https://www.gob.mx/salud/documentos/datos-abiertos-152127
Diccionarios de datos SSA https://www.gob.mx/salud/documentos/datos-abiertos-152127 y Catálogo Único de Claves de Áreas Geoestadísticas Estatales, Municipales y Localidades de INEGI https://www.inegi.org.mx/app/ageeml/
## Warning: package 'zoo' was built under R version 4.0.3
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sp
## rgdal: version: 1.5-15, (SVN revision 1045)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/vshal/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/vshal/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
ruta_estados <- '../datos_shp/Estados.shp'
estados <- rgdal::readOGR(ruta_estados)
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum International_Terrestrial_Reference_Frame_1992 in CRS
## definition: +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000
## +y_0=0 +ellps=GRS80 +units=m +no_defs
## OGR data source with driver: ESRI Shapefile
## Source: "D:\GD\Projects_actual\geo_COVID-19\datos_shp\Estados.shp", layer: "Estados"
## with 32 features
## It has 2 fields
ruta_municipios <- '../datos_shp/Municipios.shp'
municipios <- rgdal::readOGR(ruta_municipios)
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum International_Terrestrial_Reference_Frame_1992 in CRS
## definition: +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000
## +y_0=0 +ellps=GRS80 +units=m +no_defs
## OGR data source with driver: ESRI Shapefile
## Source: "D:\GD\Projects_actual\geo_COVID-19\datos_shp\Municipios.shp", layer: "Municipios"
## with 2456 features
## It has 4 fields
municipios@data$Clave_Mun_Ent_Texto <- as.numeric(as.character(municipios@data$CVE_MUN)) +
1000 * as.numeric(as.character(municipios@data$CVE_ENT))
municipios@data$Clave_Mun_Ent_Texto <- sprintf("%05d", municipios@data$Clave_Mun_Ent_Texto)
#str(municipios@data)
metadatos_m <- datos_municipio[,1:12]
#str(metadatos_m)
dia_m <- datos_municipio[,c(def_fechas)]
dia_antes_m <- cbind(datos_municipio[,c(def_fechas[-1])],data.frame(inicio = rep(0,dim(datos_municipio)[1])))
incremento_m <- dia_m - dia_antes_m
incremento_m_a100k <- round(100000 * incremento_m / datos_municipio$Pob_Municipio, 2)
incremento_m[is.na(incremento_m)] <- 0
incremento_m_a100k[is.na(incremento_m_a100k)] <- 0
#dim(incremento_m_a100k)
#head(incremento_m_a100k)
#print(casos_fechas)
#print(casos_fechas_ok)
semanas_lag <- c(seq(from = 0, to = n, by = 7))
def_breaks <- c(0,0.5,1,1.5,2,3,4,5,10,1000000)
def_breaks_labels <- c("0-0.5","0.5-1","1-1.5","1.5-2","2-3","3-4","4-5","5-10",">10")
def_raw_breaks <- c(0,1,2,3,4,5,10,15,20,1000000)
def_raw_breaks_labels <- c("1","2","3","4","5","6-10","11-15","16-20",">20")
def_indice_breaks <- c(-1,-0.75,-0.5,-0.25,-0.1,0.1,0.25,0.5,0.75,1)
def_indice_breaks_labels <- c("menos que esperado (-1)","-0.75","-0.5","-0.25","esperado (0)","0.25","0.5","0.75","mas que esperado (1)")
#for (i in 1:3) {
for (i in 1:46) {
selection_cols <- paste0("def",fecha_cadena(as.Date(fecha_inicio) + seq(semanas_lag[i]+1, semanas_lag[i]+8)))
#print(selection_cols)
suma_def_semanal <- apply(incremento_m[,selection_cols], 1, sum)
total_def_semanal <- sum(suma_def_semanal, na.rm = TRUE)
total_poblacion <- sum(datos_municipio$Pob_Municipio, na.rm = TRUE)
def_esperados_semanal <- total_def_semanal * datos_municipio$Pob_Municipio / total_poblacion
## indice (similar a NDVI)
def_indice_semanal <- (suma_def_semanal - def_esperados_semanal) / (suma_def_semanal + def_esperados_semanal)
suma_a100k_semanal <- apply(incremento_m_a100k[,selection_cols], 1, sum)
datos_semanales <- data.frame(
Clave_Mun_Ent_Texto = as.character(metadatos_m[,"Clave_Mun_Ent_Texto"]),
Municipio = metadatos_m[,"Nom_Mun"],
Latitud = metadatos_m[,"Lat_Decimal"],
Longitud = metadatos_m[,"Lon_Decimal"],
Def_100k = suma_a100k_semanal,
Def_raw = suma_def_semanal,
Def_esperado = def_esperados_semanal,
Def_indice = def_indice_semanal
)
#print(str(datos_semanales))
#print(head(datos_semanales, n = 30L))
mis_municipios <- municipios
mis_municipios <- merge(municipios, datos_semanales, by = "Clave_Mun_Ent_Texto", all.x = TRUE)
mis_municipios@data$Class_100k <- as.numeric(cut(mis_municipios@data$Def_100k,
breaks = def_breaks))
mis_municipios@data$Color_100k <- brewer.pal(n = 9, name = 'YlOrRd')[mis_municipios@data$Class_100k]
mis_municipios@data$Class_raw <- as.numeric(cut(mis_municipios@data$Def_raw,
breaks = def_raw_breaks))
mis_municipios@data$Color_raw <- brewer.pal(n = 9, name = 'PuRd')[mis_municipios@data$Class_raw]
mis_municipios@data$Class_indice <- as.numeric(cut(mis_municipios@data$Def_indice,
breaks = def_indice_breaks))
mis_municipios@data$Color_indice <- brewer.pal(n = 9, name = 'RdYlBu')[10 - mis_municipios@data$Class_indice]
#print(head(mis_municipios@data, n = 25L))
#str(mis_municipios@data)
#str(mis_municipios)
## grafica defunciones por 100 mil
plot(mis_municipios,
col = mis_municipios$Color_100k,
axes = TRUE,
border = NA,
#col = "transparent",
main = paste0("Defunciones semanales por 100 mil habitantes ",
as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
as.Date(fecha_inicio) + semanas_lag[i] + 8))
plot(estados, add = TRUE)
legend("topright",
c("Número de defunciones por COVID-19 con base en los reportes SSA de México"))
legend("bottomleft",
as.character(def_breaks_labels),
fill = brewer.pal(n = 9, name = 'YlOrRd')[1:9],
title = "Por 100 mil habitantes"
)
## grafica defunciones
plot(mis_municipios,
col = mis_municipios$Color_raw,
axes = TRUE,
border = NA,
#col = "transparent",
main = paste0("Defunciones semanales ",
as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
as.Date(fecha_inicio) + semanas_lag[i] + 8))
plot(estados, add = TRUE)
legend("topright",
c("Número de defunciones por COVID-19 con base en los reportes SSA de México"))
legend("bottomleft",
as.character(def_raw_breaks_labels),
fill = brewer.pal(n = 9, name = 'PuRd')[1:9],
title = "Defunciones semanales"
)
## grafica defuunciones indice
plot(mis_municipios,
col = mis_municipios$Color_indice,
axes = TRUE,
border = NA,
#col = "transparent",
main = paste0("Índice de desviación del número de defunciones semanales esperado ",
as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
as.Date(fecha_inicio) + semanas_lag[i] + 8))
plot(estados, add = TRUE)
legend("topright",
c("Áreas en blanco reprezentan zonas sin defunciones registrados en el periodo",
"Número de defunciones por COVID-19 con base en los reportes SSA de México"))
legend("bottomleft",
as.character(def_indice_breaks_labels),
fill = brewer.pal(n = 9, name = 'RdYlBu')[9:1],
title = "Indice de defunciones"
)
## verification plot wit coordinates only
#plot(x = mis_municipios@data$Longitud,
# y = mis_municipios@data$Latitud,
# col = mis_municipios@data$Color_raw,
# type = "p", pch = 19)
}
knitr::asis_output(htmltools::htmlPreserve(cc_html))
.